Bayesian GLM Part3

1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # for the easystats ecosystem
library(INLA) # for approximate Bayes
library(INLAutils) # for additional INLA outputs
library(patchwork) # for multiple plots
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")

2 Scenario

Here is a modified example from Peake and Quinn (1993). Peake and Quinn (1993) investigated the relationship between the number of individuals of invertebrates living in amongst clumps of mussels on a rocky intertidal shore and the area of those mussel clumps.

Figure 1: Mussels

Table 1: Format of peakquinn.csv data files

AREA	INDIV
516.00	18
469.06	60
462.25	57
938.60	100
1357.15	48
...	...

Table 2: Description of the variables in the peakquinn data file

AREA	Area of mussel clump mm2 - Predictor variable
INDIV	Number of individuals found within clump - Response variable

The aim of the analysis is to investigate the relationship between mussel clump area and the number of non-mussel invertebrate individuals supported in the mussel clump.

3 Read in the data

peake <- read_csv("../data/peakquinn.csv", trim_ws = TRUE)

Rows: 25 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): AREA, INDIV

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse

glimpse(peake)

Rows: 25
Columns: 2
$ AREA  <dbl> 516.00, 469.06, 462.25, 938.60, 1357.15, 1773.66, 1686.01, 1786.…
$ INDIV <dbl> 18, 60, 57, 100, 48, 118, 148, 214, 225, 283, 380, 278, 338, 274…

head

## Explore the first 6 rows of the data
head(peake)

# A tibble: 6 × 2
   AREA INDIV
  <dbl> <dbl>
1  516     18
2  469.    60
3  462.    57
4  939.   100
5 1357.    48
6 1774.   118

str

str(peake)

spc_tbl_ [25 × 2] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ AREA : num [1:25] 516 469 462 939 1357 ...
 $ INDIV: num [1:25] 18 60 57 100 48 118 148 214 225 283 ...
 - attr(*, "spec")=
  .. cols(
  ..   AREA = col_double(),
  ..   INDIV = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 

Easystats (datawizard)

peake |> datawizard::data_codebook()

peake (25 rows and 2 variables, 2 shown)

ID | Name  | Type    | Missings |          Values |  N
---+-------+---------+----------+-----------------+---
1  | AREA  | numeric | 0 (0.0%) | [462.25, 27144] | 25
---+-------+---------+----------+-----------------+---
2  | INDIV | numeric | 0 (0.0%) |      [18, 1402] | 25
------------------------------------------------------

Skim (modelsummary)

peake |> modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max	Histogram
AREA	25	0	7802.0	8266.8	462.2	4451.7	27144.0
INDIV	25	0	446.9	378.5	18.0	338.0	1402.0

4 Exploratory data analysis

When exploring these data as part of a frequentist analysis, exploratory data analysis revealed that the both the response (counds of individuals) and predictor (mussel clump area) were skewed and the relationship between raw counds and mussel clump area was not linear. Furthermore, there was strong evidence of a relationship between mean and variance. Normalising both reponse and predictor addressed these issues. However, rather than log transform the response, it was considered more appropriate to model against a distribution that used a logarithmic link function.

The individual observations here (\(y_i\)) are the observed number of (non mussel individuals found in mussel clump \(i\). As a count, these might be expected to follow a Poisson (or perhaps negative binomial) distribution. In the case of a negative binomial, the observed count for any given mussel clump area are expected to be drawn from a negative binomial distribution with a mean of \(\lambda_i\). All the negative binomial distributions are expected to share the same degree of dispersion (\(\theta\)) - that is, the degree to which the inhabitants of mussell clumps aggregate together (or any other reason for overdispersion) is independent of mussel clump area and can be estimated as a constant across all populations.

The natural log of the expected values (\(\lambda_i\)) is modelled against a linear predictor that includes an intercept (\(\beta_0\)) and a slope (\(\beta_1\)) associated with the natural log of mussel area.

The priors on \(\beta_0\) and \(\beta_1\) should be on the natura log scale (since this will be the scale of the parameters). As starting points, we will consider the following priors:

• \(\beta_0\): Normal prior centered at 0 with a variance of 5
• \(\beta_1\): Normal prior centered at 0 with a variance of 2
• \(\theta\): Exponential prior with rate 1

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{NB}(\lambda_i, \theta)\\ ln(\lambda_i) &= \beta_0 + \beta_1 ln(x_i)\\ \beta_0 & \sim\mathcal{N}(6,1.5)\\ \beta_1 & \sim\mathcal{N}(0,1)\\ \theta &\sim{} \mathcal{Exp}(1)\\ OR\\ \theta &\sim{} \mathcal{\Gamma}(2,1)\\ \end{align} \]

Scatterplot

ggplot(peake, aes(y = INDIV, x = AREA)) +
  geom_point() +
  geom_smooth()

Conclusions:

• there is some evidence of non-homogeneity of variance (observations are more tightly clustered around the trendline for small values of AREA and the spread increases throughout the trend).
• there is some evidence of non-linearity as evidenced by the substantial change in trajectory of the loess smoother.
• nevertheless, there is also evidence of non-normality in both the y-axis and x-axis (there are more points at low values of both y and x). Deviations from normality can make it difficult to assess other assumptions. Often the cause of homogeneity and linearity is non-homogeneity.

Boxplots

ggplot(peake, aes(y = INDIV)) +
  geom_boxplot()

ggplot(peake, aes(y = AREA)) +
  geom_boxplot()

Conclusions:

• indeed both the response (INDIV) and predictor (AREA) appear to be skewed. It is not surprising that the response is skewed as it represents counts of the number of individuals. We might expect that this should follow a Poisson distribution. Poisson distributions must be bounded by 0 at the lower end and hence as the expected value (=mean and location) of a Poisson distribution approaches 0, the distribution will become more and more asymmetrical. A the same time, the distribution would be expected to become narrower (have a smaller variance) - since in a Poisson distribution the mean and variance are equal.

Along with modelling these data against a Poisson distribution (with a natural log link), we should probably attempt to normalise the predictor variable (AREA). Whilst linear models don’t assume that the predictor variable follows a normal distribution (it is typically assumed to be a uniform distribution), it is assumed to be symmetrical. In this case, a natural log transformation might help achieve this. At the same time, it might also help homogeneity of variance and linearity.

We can mimic the action of using a log link and log-transformed predictor by presenting the scatterplot on log-transformed axes.

Scatterplot

ggplot(peake, aes(y = INDIV, x = AREA)) +
  geom_point() +
  geom_smooth() +
  scale_y_log10() +
  scale_x_log10()

Conclusions:

• there is no obvious violations of the assumptions now.

5 Fit the model

frequentist

summary(glm(INDIV ~ log(AREA), data = peake, family = poisson()))

Call:
glm(formula = INDIV ~ log(AREA), family = poisson(), data = peake)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 0.013528   0.091147   0.148    0.882    
log(AREA)   0.691628   0.009866  70.100   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 7494.5  on 24  degrees of freedom
Residual deviance: 1547.8  on 23  degrees of freedom
AIC: 1738.9

Number of Fisher Scoring iterations: 4

summary(MASS::glm.nb(INDIV ~ log(AREA), data = peake))

Call:
MASS::glm.nb(formula = INDIV ~ log(AREA), data = peake, init.theta = 7.369263003, 
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.16452    0.53387  -2.181   0.0292 *  
log(AREA)    0.82469    0.06295  13.102   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(7.3693) family taken to be 1)

    Null deviance: 161.076  on 24  degrees of freedom
Residual deviance:  25.941  on 23  degrees of freedom
AIC: 312.16

Number of Fisher Scoring iterations: 1

              Theta:  7.37 
          Std. Err.:  2.16 

 2 x log-likelihood:  -306.16 

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

peake.rstanarm <- stan_glm(INDIV ~ log(AREA),
  data = peake,
  family = poisson(),
  iter = 5000, warmup = 1000,
  chains = 3, thin = 5, refresh = 0
)

prior_summary(peake.rstanarm)

Priors for model 'peake.rstanarm' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = 0, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 0, scale = 2)
------
See help('prior_summary.stanreg') for more details

This tells us:

• for the intercept, it is using a normal prior with a mean of 0 and a standard deviation of 2.5. The are the defaults for all non-Gaussian intercepts.

• for the coefficients (in this case, just the slope), the default prior is a normal prior centered around 0 with a standard deviation of 2.5. For Poisson models, this is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the (in this case log) predictor (then rounded).

2.5 / sd(log(peake$AREA))

[1] 2.025039

• there is no auxillary parameter and thus there is no auxillary prior

Assessing priors

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refittin the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.

peake.rstanarm1 <- update(peake.rstanarm, prior_PD = TRUE)

ggemmeans(peake.rstanarm1, ~AREA) |> plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

ggemmeans(peake.rstanarm1, ~AREA) |> plot(jitter = FALSE, show_data = TRUE) + scale_y_log10()

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Conclusions:

• we see that although the range of predictions are very wide and the slope could range from strongly negative to strongly positive, the choice to have the intercept prior have a mean of 0 results in most of the prior-only posterior draws being dramatically lower than the observed values - indeed some observed values are above the range of the prior-posterior predictions.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

• \(\beta_0\): normal centred at 6 with a standard deviation of 6
  • mean of 6: because (mean(log(peake$INDIV)))
  • sd of 2.8: because (2.5 * sd(log(peake$INDIV)))
• \(\beta_1\): normal centred at 0 with a standard deviation of 1
  • sd of 2.3: because (sd(log(peake$INDIV))/sd(log(peake$AREA)))

I will also overlay the raw data for comparison.

peake.rstanarm2 <- stan_glm(INDIV ~ log(AREA),
  data = peake,
  family = poisson(),
  prior_intercept = normal(6, 2.8, autoscale = FALSE),
  prior = normal(0, 2.3, autoscale = FALSE),
  prior_PD = TRUE,
  iter = 5000, warmup = 1000,
  chains = 3, thin = 5, refresh = 0
)

ggemmeans(peake.rstanarm2, ~AREA) |>
  plot(show_data = TRUE) +
  scale_y_log10()

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Now lets refit, conditioning on the data.

peake.rstanarm3 <- update(peake.rstanarm2, prior_PD = FALSE)

Plotting prior and posterior

posterior_vs_prior(peake.rstanarm3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Drawing from prior...

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

ggemmeans(peake.rstanarm3, ~AREA) |> plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

form <- bf(INDIV ~ log(AREA), family = poisson())
peake.brm <- brm(form,
  data = peake,
  iter = 5000, warmup = 1000,
  chains = 3, cores = 3,
  thin = 5, refresh = 0,
  backend = "rstan"
)

Compiling Stan program...

Start sampling

prior_summary(peake.brm)

                  prior     class    coef group resp dpar nlpar lb ub       source
                 (flat)         b                                          default
                 (flat)         b logAREA                             (vectorized)
 student_t(3, 5.8, 2.5) Intercept                                          default

This tells us:

• for the intercept, it is using a student t (flatter normal) prior with a mean of 5.8 and a standard deviation of 2.5. The mean of this prior is based on the median of the log-transformed response and the standard deviation of 2.5 is the default minimum.

[1] 5.823046

mad(log(peake$INDIV))

[1] 1.025465

• for the beta coefficients (in this case, just the slope), the default prior is a inproper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.

• there is no sigma parameter in a binomial model and therefore there are no additional priors.

Assessing priors

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior='only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.

We will adopt a similar approach to that of rstanarm - that is 2.5 * sd(log(peake$INDIV))/sd(log(peake$AREA)) (approx. 2.3).

form <- bf(INDIV ~ log(AREA), family = poisson())
peake.brm1 <- brm(form,
  data = peake,
  prior = c(
    prior(normal(0, 2.3), class = "b")
  ),
  sample_prior = "only",
  iter = 5000, warmup = 1000,
  chains = 3, cores = 3,
  thin = 5, refresh = 0,
  backend = "rstan"
)

Compiling Stan program...

Start sampling

ggemmeans(peake.brm1, ~AREA) |> plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

ggemmeans(peake.brm1, ~AREA) |> plot(show_data = TRUE) + scale_y_log10()

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

conditional_effects(peake.brm1) |>
  plot(points = TRUE) |>
  _[[1]] + scale_y_log10()

plot(conditional_effects(peake.brm1), points = TRUE, plot = FALSE)[[1]] + scale_y_log10()

Conclusions:

• we see that the range of predictions is faily wide and the slope could range from strongly negative to strongly positive.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

• \(\beta_0\): normal centred at 6 with a standard deviation of 6
  • mean of 6: because (mean(log(peake$INDIV)))
  • sd of 1.5: because (sd(log(peake$INDIV)))
• \(\beta_1\): normal centred at 0 with a standard deviation of 1
  • sd of 1: because (sd(log(peake$INDIV))/sd(log(peake$AREA)))

I will also overlay the raw data for comparison.

form <- bf(INDIV ~ log(AREA), family = poisson())
priors <- prior(normal(6, 1.5), class = "Intercept") +
  prior(normal(0, 1), class = "b")
peake.brm2 <- brm(form,
  data = peake,
  prior = priors,
  sample_prior = "only",
  iter = 5000, warmup = 1000,
  chains = 3, cores = 3,
  thin = 5, refresh = 0,
  backend = "rstan"
)

Compiling Stan program...

Start sampling

ggemmeans(peake.brm2, ~AREA) |>
  plot(show_data = TRUE) +
  scale_y_log10()

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

peake.brm3 <- update(peake.brm2, sample_prior = TRUE, refresh = 0)

The desired updates require recompiling the model

Compiling Stan program...

Start sampling

Plotting prior and posterior

peake.brm3 |> get_variables()

 [1] "b_Intercept"     "b_logAREA"       "Intercept"       "prior_Intercept"
 [5] "prior_b"         "lprior"          "lp__"            "accept_stat__"  
 [9] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
[13] "energy__"       

peake.brm3 |>
  hypothesis("logAREA=0") |>
  plot()

peake.brm3 |> SUYR_prior_and_posterior()

Exploring the stan code

standata(peake.brm3)

$N
[1] 25

$Y
 [1]   18   60   57  100   48  118  148  214  225  283  380  278  338  274  569
[16]  509  682  600  978  363 1402  675 1159 1062  632

$K
[1] 2

$Kc
[1] 1

$X
   Intercept   logAREA
1          1  6.246107
2          1  6.150731
3          1  6.136106
4          1  6.844389
5          1  7.213142
6          1  7.480800
7          1  7.430120
8          1  7.487896
9          1  8.035949
10         1  8.289067
11         1  8.394989
12         1  8.401037
13         1  8.513765
14         1  8.400853
15         1  8.610818
16         1  8.919481
17         1  8.873303
18         1  9.121503
19         1  9.223560
20         1  9.136445
21         1  9.527275
22         1  9.918414
23         1 10.115073
24         1 10.208912
25         1 10.170373
attr(,"assign")
[1] 0 1

$prior_only
[1] 0

attr(,"class")
[1] "standata" "list"    

stancode(peake.brm3)

// generated with brms 2.22.0
functions {
}
data {
  int<lower=1> N;  // total number of observations
  array[N] int Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int<lower=1> Kc;  // number of population-level effects after centering
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
}
parameters {
  vector[Kc] b;  // regression coefficients
  real Intercept;  // temporary intercept for centered predictors
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += normal_lpdf(b | 0, 1);
  lprior += normal_lpdf(Intercept | 6, 1.5);
}
model {
  // likelihood including constants
  if (!prior_only) {
    target += poisson_log_glm_lpmf(Y | Xc, Intercept, b);
  }
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
  // additionally sample draws from priors
  real prior_b = normal_rng(0,1);
  real prior_Intercept = normal_rng(6,1.5);
}

6 MCMC sampling diagnostics

rstanarm

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

See list of available diagnostics by name

available_mcmc()

bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

• mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

plot(peake.rstanarm3, plotfun = "mcmc_trace")

The chains appear well mixed and very similar

• acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

plot(peake.rstanarm3, "acf_bar")

There is no evidence of autocorrelation in the MCMC samples

• Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

plot(peake.rstanarm3, "rhat_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

plot(peake.rstanarm3, "neff_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

More diagnostics

plot(peake.rstanarm3, "combo")

plot(peake.rstanarm3, "violin")

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(peake.rstanarm3)

The chains appear well mixed and very similar

• stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(peake.rstanarm3)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(peake.rstanarm3)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

stan_ess(peake.rstanarm3)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

stan_dens(peake.rstanarm3, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

peake.ggs <- ggs(peake.rstanarm3)
ggs_traceplot(peake.ggs)

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(peake.ggs)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(peake.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

ggs_effective(peake.ggs)

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
  always returns an ungrouped data frame and adjust accordingly.
ℹ The deprecated feature was likely used in the ggmcmc package.
  Please report the issue at <https://github.com/xfim/ggmcmc/issues/>.

Ratios all very high.

More diagnostics

ggs_crosscorrelation(peake.ggs)

ggs_grb(peake.ggs)

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

See list of available diagnostics by name

available_mcmc()

bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

• trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

mcmc_plot(peake.brm3, type = "trace")

No divergences to plot.

The chains appear well mixed and very similar

• acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

mcmc_plot(peake.brm3, type = "acf_bar")

There is no evidence of autocorrelation in the MCMC samples

• rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mcmc_plot(peake.brm3, type = "rhat_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

mcmc_plot(peake.brm3, type = "neff_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

More diagnostics

mcmc_plot(peake.brm3, type = "combo")

mcmc_plot(peake.brm3, type = "violin")

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(peake.brm3$fit)

The chains appear well mixed and very similar

• stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(peake.brm3$fit)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(peake.brm3$fit)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

stan_ess(peake.brm3$fit)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

stan_dens(peake.brm3$fit, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

peake.ggs <- ggs(peake.brm3)
ggs_traceplot(peake.ggs)

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(peake.ggs)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(peake.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

ggs_effective(peake.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(peake.ggs)

ggs_grb(peake.ggs)

7 Model validation

rstanarm

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

See list of available diagnostics by name

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(peake.rstanarm3, plotfun = "dens_overlay")

The model draws appear deviate from the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

pp_check(peake.rstanarm3, plotfun = "error_scatter_avg")

The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.

• error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. Again, this is not interpretable for binary data.

pp_check(peake.rstanarm3, x = peake$AREA, plotfun = "error_scatter_avg_vs_x")

• intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(peake.rstanarm3, x = peake$AREA, plotfun = "intervals")

The modelled preditions seem to underestimate the uncertainty with increasing mussel clump area.

• ribbon: this is just an alternative way of expressing the above plot.

pp_check(peake.rstanarm3, x = peake$AREA, plotfun = "ribbon")

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(peake.rstanarm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(peake.rstanarm3, ndraws = 250, summary = FALSE)
peake.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = peake$INDIV,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(peake.resids)

peake.resids |> testDispersion()

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 97.324, p-value < 2.2e-16
alternative hypothesis: two.sided

Conclusions:

• the simulated residuals suggest a general lack of fit due to overdispersion and outliers
• perhaps we should explore a negative binomial model

brms

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

See list of available diagnostics by name

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(peake.brm3, type = "dens_overlay", ndraws = 100)

The model draws appear deviate from the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

pp_check(peake.brm3, type = "error_scatter_avg")

Using all posterior draws for ppc type 'error_scatter_avg' by default.

The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.

• error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. Again, this is not interpretable for binary data.

pp_check(peake.brm3, x = "AREA", type = "error_scatter_avg_vs_x")

Using all posterior draws for ppc type 'error_scatter_avg_vs_x' by default.

• intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(peake.brm3, x = "AREA", type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

The modelled preditions seem to underestimate the uncertainty with increasing mussel clump area.

• ribbon: this is just an alternative way of expressing the above plot.

pp_check(peake.brm3, x = "AREA", type = "ribbon")

Using all posterior draws for ppc type 'ribbon' by default.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(peake.brm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(peake.brm3, ndraws = 250, summary = FALSE)
peake.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = peake$INDIV,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(peake.resids)

peake.resids |> testDispersion()

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 98.867, p-value < 2.2e-16
alternative hypothesis: two.sided

peake.resids <- make_brms_dharma_res(peake.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(peake.resids)) +
  wrap_elements(~ plotResiduals(peake.resids, form = factor(rep(1, nrow(peake))))) +
  wrap_elements(~ plotResiduals(peake.resids, quantreg = FALSE)) +
  wrap_elements(~ testDispersion(peake.resids))

Conclusions:

• the simulated residuals suggest a general lack of fit due to overdispersion and outliers
• perhaps we should explore a negative binomial model

8 Fit Negative Binomial rstanarm

rstanarm

Fit the model

peake.rstanarm4 <- stan_glm(INDIV ~ log(AREA),
  data = peake,
  family = neg_binomial_2(),
  prior_intercept = normal(6, 2.8, autoscale = FALSE),
  prior = normal(0, 2.3, autoscale = FALSE),
  prior_aux = rstanarm::exponential(rate = 1, autoscale = FALSE),
  iter = 5000, warmup = 1000,
  chains = 3, thin = 5, refresh = 0
)

Plotting prior and posterior

posterior_vs_prior(peake.rstanarm4,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Drawing from prior...

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

## ggemmeans(peake.rstanarm4,  ~AREA) |> plot(show_data=TRUE)
ggpredict(peake.rstanarm4, ~AREA) |> plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

MCMC chain and posterior checks

plot(peake.rstanarm4, plotfun = "mcmc_trace")

plot(peake.rstanarm4, "acf_bar")

plot(peake.rstanarm4, "rhat_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot(peake.rstanarm4, "neff_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(peake.rstanarm4, x = peake$AREA, plotfun = "dens_overlay")

Warning: The following arguments were unrecognized and ignored: x

pp_check(peake.rstanarm4, x = peake$AREA, plotfun = "error_scatter_avg_vs_x")

pp_check(peake.rstanarm4, x = peake$AREA, plotfun = "intervals")

DHARMa residuals

preds <- posterior_predict(peake.rstanarm4, ndraws = 250, summary = FALSE)
peake.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = peake$INDIV,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(peake.resids)

LOOIC

Different models can be compared via information criterion. LOO (Leave One Out) IC is similar to AIC except that AIC does not consider priors and assumes that the posterior likelihood is multivariate normal. LOO AIC does not and integrates over all uncertainty.

• ELPD: is the theoretical expected log pointwise predictive density for a new dataset and can be estimated via leave-one-out cross-validation.
• elpd_loo: is the Bayesian LOO estimate of ELPD
• SE of elpd_loo: is the Monte Carlo standard error is the estimate for the computational accuracy of MCMC and importance sampling used to compute elpd_loo
• p_loo: is the difference between elpd_loo and the estimated from non-cross validated ELPD and therefore is a measure of the relative extra difficulty of predicting new data compared to the observed data. It can also be a meaure of the effective number of parameters:
  • in a well behaved model, p_loo will be less than the number of estimated parameters and the total sample size.
  • in a misspecified model, p_loo will be greater than the number of estimated parameters and the total sample size.
• Pareto k: is a relative measure of the influence of each observation as well as the accuracy with which the associated leave one out calculations can be estimated.
  • when k < 0.5: the corresponding components can be accurately estimated
  • when 0.5 < k < 0.7: the accuracy is lower but still acceptable
  • when k>0.7: the accuracy is too low and elpd_loo is unreliable. This can also suggest a misspecified model.

More information about the above can be gleaned from [https://mc-stan.org/loo/reference/loo-glossary.html].

Start with the Poisson model

(peake.rstanarm3.loo <- loo(peake.rstanarm3))

Warning: Found 5 observation(s) with a pareto_k > 0.7. We recommend calling 'loo' again with argument 'k_threshold = 0.7' in order to calculate the ELPD without the assumption that these observations are negligible. This will refit the model 5 times to compute the ELPDs for the problematic observations directly.

Computed from 2400 by 25 log-likelihood matrix.

         Estimate    SE
elpd_loo   -930.3 298.1
p_loo       113.5  43.8
looic      1860.7 596.2
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     20    80.0%   69      
   (0.7, 1]   (bad)       1     4.0%   <NA>    
   (1, Inf)   (very bad)  4    16.0%   <NA>    
See help('pareto-k-diagnostic') for details.

Conclusions:

• p_loo is substantially higher than the number of estimated parameters
• there are some Pareto k values identified as bad or even very bad

Now for the Negative Binomial model

(peake.rstanarm4.loo <- loo(peake.rstanarm4))

Computed from 2400 by 25 log-likelihood matrix.

         Estimate   SE
elpd_loo   -156.7  5.3
p_loo         1.9  0.5
looic       313.5 10.6
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

Conclusions:

• p_loo is lower than the number of estimated parameters
• there are no bad Pareto k values

We can also compare the models.

The difference in Expected Log Pointwise Probability Density (elpd_dff) computes the difference in elpd_loo of two models (expressed relative to the higher model). The difference in elpd_diff will be negative if the expected out-of-sample predictive accuracy of the first model is higher. If the difference is be positive then the second model is preferred.

Note, the above assumes that both models are valid. In the current case, we know that the Poisson model was overdispersed and thus it should not be a genuine candidate. Nevertheless, we will compare the Poisson to the Negative Binomial for illustrative purposes.

loo_compare(peake.rstanarm3.loo, peake.rstanarm4.loo)

                elpd_diff se_diff
peake.rstanarm4    0.0       0.0 
peake.rstanarm3 -773.6     294.5 

Conclusions:

• the elpd_diff is negative, indicating that the first model (Negative Binomial) is better.

brms

Fit the model

form <- bf(INDIV ~ log(AREA), family = negbinomial())
priors <-
  prior(normal(6, 1.5), class = "Intercept") +
  prior(normal(0, 1), class = "b") +
  ## prior(gamma(0.01, 0.01), class='shape')
  prior(gamma(2, 1), class = "shape")
peake.brm4 <- brm(form,
  data = peake,
  prior = priors,
  sample_prior = TRUE,
  iter = 5000, warmup = 1000,
  chains = 3, cores = 3,
  thin = 5, refresh = 0,
  backend = "rstan"
)

Compiling Stan program...

Start sampling

Plotting prior and posterior

peake.brm4 |> SUYR_prior_and_posterior()

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

## ggemmeans(peake.brm4,  ~AREA) |> plot(show_data=TRUE)
ggpredict(peake.brm4, ~AREA) |> plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

MCMC chain and posterior checks

peake.brm4$fit |> stan_trace()

peake.brm4$fit |> stan_ac()

peake.brm4$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

peake.brm4$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(peake.brm4, x = "AREA", type = "dens_overlay", ndraws = 100)

Warning: The following arguments were unrecognized and ignored: x

pp_check(peake.brm4, x = "AREA", type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

DHARMa residuals

preds <- posterior_predict(peake.brm4, ndraws = 250, summary = FALSE)
peake.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = peake$INDIV,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(peake.resids)

peake.resids <- make_brms_dharma_res(peake.brm4, integerResponse = FALSE)
wrap_elements(~ testUniformity(peake.resids)) +
  wrap_elements(~ plotResiduals(peake.resids, form = factor(rep(1, nrow(peake))))) +
  wrap_elements(~ plotResiduals(peake.resids, quantreg = FALSE)) +
  wrap_elements(~ testDispersion(peake.resids))

LOOIC

Different models can be compared via information criterion. LOO (Leave One Out) IC is similar to AIC except that AIC does not consider priors and assumes that the posterior likelihood is multivariate normal. LOO AIC does not and integrates over all uncertainty.

• ELPD: is the theoretical expected log pointwise predictive density for a new dataset and can be estimated via leave-one-out cross-validation.
• elpd_loo: is the Bayesian LOO estimate of ELPD
• SE of elpd_loo: is the Monte Carlo standard error is the estimate for the computational accuracy of MCMC and importance sampling used to compute elpd_loo
• p_loo: is the difference between elpd_loo and the estimated from non-cross validated ELPD and therefore is a measure of the relative extra difficulty of predicting new data compared to the observed data. It can also be a meaure of the effective number of parameters:
  • in a well behaved model, p_loo will be less than the number of estimated parameters and the total sample size.
  • in a misspecified model, p_loo will be greater than the number of estimated parameters and the total sample size.
• Pareto k: is a relative measure of the influence of each observation as well as the accuracy with which the associated leave one out calculations can be estimated.
  • when k < 0.5: the corresponding components can be accurately estimated
  • when 0.5 < k < 0.7: the accuracy is lower but still acceptable
  • when k>0.7: the accuracy is too low and elpd_loo is unreliable. This can also suggest a misspecified model.

More information about the above can be gleaned from [https://mc-stan.org/loo/reference/loo-glossary.html].

Start with the Poisson model

(peake.brm3.loo <- loo(peake.brm3))

Warning: Found 8 observations with a pareto_k > 0.7 in model 'peake.brm3'. We
recommend to set 'moment_match = TRUE' in order to perform moment matching for
problematic observations.

Computed from 2400 by 25 log-likelihood matrix.

         Estimate    SE
elpd_loo   -924.5 295.1
p_loo       108.4  41.1
looic      1849.0 590.3
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     17    68.0%   75      
   (0.7, 1]   (bad)       4    16.0%   <NA>    
   (1, Inf)   (very bad)  4    16.0%   <NA>    
See help('pareto-k-diagnostic') for details.

Conclusions:

• p_loo is substantially higher than the number of estimated parameters
• there are some Pareto k values identified as bad or even very bad

Now for the Negative Binomial model

(peake.brm4.loo <- loo(peake.brm4))

Computed from 2400 by 25 log-likelihood matrix.

         Estimate   SE
elpd_loo   -156.4  5.4
p_loo         2.0  0.5
looic       312.8 10.8
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.0]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

Conclusions:

• p_loo is lower than the number of estimated parameters
• there are no bad Pareto k values

We can also compare the models.

The difference in Expected Log Pointwise Probability Density (elpd_dff) computes the difference in elpd_loo of two models (expressed relative to the higher model). The difference in elpd_diff will be negative if the expected out-of-sample predictive accuracy of the first model is higher. If the difference is be positive then the second model is preferred.

Note, the above assumes that both models are valid. In the current case, we know that the Poisson model was overdispersed and thus it should not be a genuine candidate. Nevertheless, we will compare the Poisson to the Negative Binomial for illustrative purposes.

loo_compare(peake.brm3.loo, peake.brm4.loo)

           elpd_diff se_diff
peake.brm4    0.0       0.0 
peake.brm3 -768.1     291.4 

Conclusions:

• the elpd_diff is negative, indicating that the first model (Negative Binomial) is better.

9 Partial effects plots

rstanarm

ggpredict

peake.rstanarm4 |>
  ggpredict() |>
  plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

ggemmeans

## Note, does not seem to backtransform...
peake.rstanarm4 |>
  ggemmeans(~AREA, type = "fixed", back.transform = TRUE) |>
  plot()

epred_draws

peake.rstanarm4 |>
  epred_draws(newdata = peake) |>
  median_hdci() |>
  ggplot(aes(x = AREA, y = .epred)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
  geom_line() +
  geom_point(data = peake, aes(y = INDIV, x = AREA))

brms

conditional_effects

peake.brm4 |>
  conditional_effects() |>
  plot(points = TRUE)

peake.brm4 |>
  conditional_effects(spaghetti = TRUE, ndraws = 200) |>
  plot(points = TRUE) + theme_classic()

NULL

ce <- peake.brm4 |> conditional_effects(spaghetti = TRUE, ndraws = 200)
plot(ce, points = TRUE, plot = FALSE)[[1]] + theme_classic()

ggpredict

peake.brm4 |>
  ggpredict() |>
  plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

ggemmeans

peake.brm4 |>
  ggemmeans(~AREA) |>
  plot(show_data = TRUE)

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

epred_draws

peake.brm4 |>
  epred_draws(newdata = peake) |>
  median_hdci() |>
  ggplot(aes(x = AREA, y = .epred)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.3) +
  geom_line() +
  geom_point(data = peake, aes(y = INDIV, x = AREA))

10 Model investigation

rstanarm

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(peake.rstanarm4)

Model Info:
 function:     stan_glm
 family:       neg_binomial_2 [log]
 formula:      INDIV ~ log(AREA)
 algorithm:    sampling
 sample:       2400 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 25
 predictors:   2

Estimates:
                        mean   sd   10%   50%   90%
(Intercept)           -1.1    0.7 -2.0  -1.2  -0.2 
log(AREA)              0.8    0.1  0.7   0.8   0.9 
reciprocal_dispersion  4.6    1.3  3.1   4.5   6.4 

Fit Diagnostics:
           mean   sd    10%   50%   90%
mean_PPD 480.4   91.0 376.9 470.8 594.2

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
                      mcse Rhat n_eff
(Intercept)           0.0  1.0  2491 
log(AREA)             0.0  1.0  2533 
reciprocal_dispersion 0.0  1.0  2003 
mean_PPD              2.0  1.0  2159 
log-posterior         0.0  1.0  2043 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusions:

• in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 19 raw observations.
• the estimated intercept (expected number of individuals for a mussel clump of log-transformed area = 0) is -1.15. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.32. When the mussel clump is 0, we expect to find 0.32 individuals in the mussel clump (which obviously does not exist if it has a log-transformed area of 0!).
  
• the estimated slope (rate at which the number of individuals changes per 1 unit change in log-transformed area), is 0.82 (mean) or 0.72 (median) with a standard deviation of 0. The 90% credibility intervals indicate that we are 90% confident that the slope is between 0.82 and 0.82 - e.g. there is a significant positive trend. When back-transformed onto the reponse scale, we see that for a one unit increase in log-transformed mussel clump area, the number of individuals increases by a factor of 2.28. This represents a 128% increase per 1 unit change in log-transformed mussel clump area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

summarise_draws (posterior)

peake.rstanarm4$stanfit |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 5 × 8
  variable               median    lower    upper  rhat length ess_bulk ess_tail
  <chr>                   <dbl>    <dbl>    <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Intercept)            -1.17    -2.51     0.333 1.00    2400    2510.    2408.
2 log(AREA)               0.824    0.645    0.980 1.000   2400    2557.    2372.
3 reciprocal_dispersi…    4.49     2.42     7.27  1.00    2400    2041.    1940.
4 mean_PPD              471.     326.     668.    1.000   2400    2278.    2328.
5 log-posterior        -161.    -164.    -160.    1.00    2400    2240.    2227.

We can also alter the CI level.

peake.rstanarm4$stanfit |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x, credMass = 0.9),
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 5 × 8
  variable               median    lower    upper  rhat length ess_bulk ess_tail
  <chr>                   <dbl>    <dbl>    <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Intercept)            -1.17    -2.33   2.04e-2 1.00    2400    2510.    2408.
2 log(AREA)               0.824    0.680  9.58e-1 1.000   2400    2557.    2372.
3 reciprocal_dispersi…    4.49     2.54   6.67e+0 1.00    2400    2041.    1940.
4 mean_PPD              471.     336.     6.14e+2 1.000   2400    2278.    2328.
5 log-posterior        -161.    -163.    -1.60e+2 1.00    2400    2240.    2227.

Arguably, it would be better to back-transform to the ratio scale

peake.rstanarm4$stanfit |>
  ## tidy_draws() |>
  ## exp() |>
  summarise_draws(
    ~ median(exp(.x)),
    ~ HDInterval::hdi(exp(.x)),
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 5 × 8
  variable `~median(exp(.x))`     lower     upper  rhat length ess_bulk ess_tail
  <chr>                 <dbl>     <dbl>     <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Interc…          3.11e-  1 4.11e-  2 1.09e+  0 1.00    2400    2510.    2408.
2 log(ARE…          2.28e+  0 1.89e+  0 2.65e+  0 1.000   2400    2557.    2372.
3 recipro…          8.95e+  1 4.89e+  0 1.06e+  3 1.00    2400    2041.    1940.
4 mean_PPD          2.87e+204 1.06e+102 3.81e+279 1.000   2400    2278.    2328.
5 log-pos…          7.63e- 71 3.73e- 75 2.21e- 70 1.00    2400    2240.    2227.

tidyMCMC

tidyMCMC(peake.rstanarm4$stanfit, estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval", rhat = TRUE, ess = TRUE)

# A tibble: 5 × 7
  term                  estimate std.error conf.low conf.high  rhat   ess
  <chr>                    <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
1 (Intercept)             -1.15     0.724    -2.51      0.333 1.000  2491
2 log(AREA)                0.824    0.0850    0.645     0.980 1.000  2533
3 reciprocal_dispersion    4.64     1.29      2.42      7.27  1.000  2003
4 mean_PPD               480.      91.0     326.      668.    0.999  2159
5 log-posterior         -162.       1.27   -164.     -160.    1.00   2043

Conclusions:

• the estimated intercept (expected number of individuals for a mussel clump of log-transformed area = 0) is -1.15. This is the median of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.32. When the mussel clump is 0, we expect to find 0.32 individuals in the mussel clump (which obviously does not exist if it has a log-transformed area of 0!).
  
• the estimated slope (rate at which the number of individuals changes per 1 unit change in log-transformed area), is 0.82 (mean) with a standard error of 0.08. The 95% credibility intervals indicate that we are 90% confident that the slope is between 0.65 and 0.98 - e.g. there is a significant positive trend. When back-transformed onto the reponse scale, we see that for a one unit increase in log-transformed mussel clump area, the number of individuals increases by a factor of 2.28. This represents a 128% increase per 1 unit change in log-transformed mussel clump area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

as_draws_df (posterior)

peake.rstanarm4$stanfit |> as_draws_df()

# A draws_df: 800 iterations, 3 chains, and 5 variables
   (Intercept) log(AREA) reciprocal_dispersion mean_PPD log-posterior
1       -0.129      0.70                   3.4      402          -162
2       -0.368      0.74                   4.0      517          -161
3       -0.969      0.82                   8.5      451          -164
4        0.057      0.68                   3.8      416          -162
5       -3.028      1.02                   2.7      337          -166
6       -1.424      0.88                   3.4      693          -162
7       -1.837      0.88                   2.8      490          -164
8       -2.082      0.91                   3.6      336          -163
9       -0.688      0.76                   3.0      416          -162
10      -1.529      0.89                   4.2      570          -163
# ... with 2390 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

## summarised
peake.rstanarm4$stanfit |>
  as_draws_df() |>
  exp() |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    ess_bulk, ess_tail
  )

# A tibble: 5 × 7
  variable                 median     lower     upper  rhat ess_bulk ess_tail
  <chr>                     <dbl>     <dbl>     <dbl> <dbl>    <dbl>    <dbl>
1 (Intercept)           3.11e-  1 4.11e-  2 1.09e+  0 1.00     2510.    2408.
2 log(AREA)             2.28e+  0 1.89e+  0 2.65e+  0 1.000    2557.    2372.
3 reciprocal_dispersion 8.95e+  1 4.89e+  0 1.06e+  3 1.000    2041.    1940.
4 mean_PPD              2.87e+204 1.06e+102 3.81e+279 1.000    2302.    2328.
5 log-posterior         7.63e- 71 3.73e- 75 2.21e- 70 1.00     2240.    2227.

## summarised on fractional scale
peake.rstanarm4$stanfit |>
  as_draws_df() |>
  dplyr::select(matches("Intercept|AREA")) |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    length,
    ess_bulk, ess_tail
  )

Warning: Dropping 'draws_df' class as required metadata was removed.

# A tibble: 2 × 8
  variable    median  lower upper  rhat length ess_bulk ess_tail
  <chr>        <dbl>  <dbl> <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Intercept) -1.17  -2.51  0.333 1.000   2400    2479.    2372.
2 log(AREA)    0.824  0.645 0.980 1.000   2400    2526.    2341.

## OR
names <- peake.rstanarm4 |>
  formula() |>
  model.matrix(peake) |>
  colnames()

peake.rstanarm4$stanfit |>
  as_draws_df() |>
  dplyr::select(any_of(names)) |>
  mutate(across(everything(), exp)) |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    length,
    ess_bulk, ess_tail
  )

Warning: Dropping 'draws_df' class as required metadata was removed.

# A tibble: 2 × 8
  variable    median  lower upper  rhat length ess_bulk ess_tail
  <chr>        <dbl>  <dbl> <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 (Intercept)  0.311 0.0411  1.09 1.000   2400    2479.    2372.
2 log(AREA)    2.28  1.89    2.65 1.000   2400    2526.    2341.

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

peake.rstanarm4 |> get_variables()

[1] "(Intercept)"           "log(AREA)"             "reciprocal_dispersion"
[4] "accept_stat__"         "stepsize__"            "treedepth__"          
[7] "n_leapfrog__"          "divergent__"           "energy__"             

peake.draw <- peake.rstanarm4 |> gather_draws(`.Intercept.*|.*AREA.*`, regex = TRUE)
peake.draw

# A tibble: 4,800 × 5
# Groups:   .variable [2]
   .chain .iteration .draw .variable    .value
    <int>      <int> <int> <chr>         <dbl>
 1      1          1     1 (Intercept) -0.129 
 2      1          2     2 (Intercept) -0.368 
 3      1          3     3 (Intercept) -0.969 
 4      1          4     4 (Intercept)  0.0566
 5      1          5     5 (Intercept) -3.03  
 6      1          6     6 (Intercept) -1.42  
 7      1          7     7 (Intercept) -1.84  
 8      1          8     8 (Intercept) -2.08  
 9      1          9     9 (Intercept) -0.688 
10      1         10    10 (Intercept) -1.53  
# ℹ 4,790 more rows

We can then summarise this

peake.draw |> median_hdci()

# A tibble: 2 × 7
  .variable   .value .lower .upper .width .point .interval
  <chr>        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 (Intercept) -1.17  -2.44   0.415   0.95 median hdci     
2 log(AREA)    0.824  0.645  0.980   0.95 median hdci     

peake.rstanarm4 |>
  gather_draws(`.Intercept.*|.*AREA.*`, regex = TRUE) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

We could alternatively express the parameters on the response scale.

peake.rstanarm4 |>
  gather_draws(`.Intercept.*|.*AREA.*`, regex = TRUE) |>
  group_by(.variable) |>
  mutate(.value = exp(.value)) |>
  median_hdci()

# A tibble: 2 × 7
  .variable   .value .lower .upper .width .point .interval
  <chr>        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 (Intercept)  0.311 0.0411   1.09   0.95 median hdci     
2 log(AREA)    2.28  1.89     2.65   0.95 median hdci     

Conclusions:

• the estimated intercept (expected number of individuals for a mussel clump of log-transformed area = 0) is -1.17. This is the median of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.31. When the mussel clump is 0, we expect to find 0.31 individuals in the mussel clump (which obviously does not exist if it has a log-transformed area of 0!).
  
• the estimated slope (rate at which the number of individuals changes per 1 unit change in log-transformed area), is 0.82 (mean) with a standard error of 0.65. The 95% credibility intervals indicate that we are 90% confident that the slope is between 0.98 and 0.95 - e.g. there is a significant positive trend. When back-transformed onto the reponse scale, we see that for a one unit increase in log-transformed mussel clump area, the number of individuals increases by a factor of 2.28. This represents a 128% increase per 1 unit change in log-transformed mussel clump area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

bayesplot

peake.rstanarm4 |> plot(plotfun = "mcmc_intervals")

tidy_draws

This is purely a graphical depiction on the posteriors.

peake.rstanarm4 |> tidy_draws()

# A tibble: 2,400 × 12
   .chain .iteration .draw `(Intercept)` `log(AREA)` reciprocal_dispersion
    <int>      <int> <int>         <dbl>       <dbl>                 <dbl>
 1      1          1     1       -0.129        0.703                  3.41
 2      1          2     2       -0.368        0.735                  3.99
 3      1          3     3       -0.969        0.816                  8.49
 4      1          4     4        0.0566       0.682                  3.84
 5      1          5     5       -3.03         1.02                   2.66
 6      1          6     6       -1.42         0.876                  3.36
 7      1          7     7       -1.84         0.881                  2.76
 8      1          8     8       -2.08         0.915                  3.55
 9      1          9     9       -0.688        0.763                  3.01
10      1         10    10       -1.53         0.892                  4.18
# ℹ 2,390 more rows
# ℹ 6 more variables: accept_stat__ <dbl>, stepsize__ <dbl>, treedepth__ <dbl>,
#   n_leapfrog__ <dbl>, divergent__ <dbl>, energy__ <dbl>

spread_draws

peake.rstanarm4 |> spread_draws(`.Intercept.*|.*AREA.*`, regex = TRUE)

# A tibble: 2,400 × 5
   .chain .iteration .draw `(Intercept)` `log(AREA)`
    <int>      <int> <int>         <dbl>       <dbl>
 1      1          1     1       -0.129        0.703
 2      1          2     2       -0.368        0.735
 3      1          3     3       -0.969        0.816
 4      1          4     4        0.0566       0.682
 5      1          5     5       -3.03         1.02 
 6      1          6     6       -1.42         0.876
 7      1          7     7       -1.84         0.881
 8      1          8     8       -2.08         0.915
 9      1          9     9       -0.688        0.763
10      1         10    10       -1.53         0.892
# ℹ 2,390 more rows

## summarised
peake.rstanarm4 |>
  spread_draws(`.Intercept.*|.*AREA.*`, regex = TRUE) |>
  summarise_draws(
    "median",
    ~ HDInterval::hdi(.x),
    "rhat",
    "ess_bulk"
  )

# A tibble: 2 × 6
  variable    median  lower upper  rhat ess_bulk
  <chr>        <dbl>  <dbl> <dbl> <dbl>    <dbl>
1 (Intercept) -1.17  -2.51  0.333 1.00     2510.
2 log(AREA)    0.824  0.645 0.980 1.000    2557.

## summarised on fractional scale
peake.rstanarm4 |>
  spread_draws(`.Intercept.*|.*AREA.*`, regex = TRUE) |>
  mutate(across(everything(), exp)) |>
  summarise_draws(
    "median",
    ~ HDInterval::hdi(.x),
    "rhat",
    "ess_bulk"
  )

# A tibble: 2 × 6
  variable    median  lower upper  rhat ess_bulk
  <chr>        <dbl>  <dbl> <dbl> <dbl>    <dbl>
1 (Intercept)  0.311 0.0411  1.09 1.00     2510.
2 log(AREA)    2.28  1.89    2.65 1.000    2557.

posterior_samples

peake.rstanarm4 |>
  posterior_samples() |>
  as_tibble()

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

# A tibble: 2,400 × 3
   `(Intercept)` `log(AREA)` reciprocal_dispersion
           <dbl>       <dbl>                 <dbl>
 1       -0.129        0.703                  3.41
 2       -0.368        0.735                  3.99
 3       -0.969        0.816                  8.49
 4        0.0566       0.682                  3.84
 5       -3.03         1.02                   2.66
 6       -1.42         0.876                  3.36
 7       -1.84         0.881                  2.76
 8       -2.08         0.915                  3.55
 9       -0.688        0.763                  3.01
10       -1.53         0.892                  4.18
# ℹ 2,390 more rows

\(R^2\)

Unfortunately, \(R^2\) calculations for models other than Gaussian and Binomial have not yet been implemented for rstanarm models yet.

## peake.rstanarm4 |> bayes_R2() |> median_hdci

brms

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(peake.brm4)

 Family: negbinomial 
  Links: mu = log; shape = identity 
Formula: INDIV ~ log(AREA) 
   Data: peake (Number of observations: 25) 
  Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
         total post-warmup draws = 2400

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -1.10      0.69    -2.47     0.28 1.00     2401     2461
logAREA       0.82      0.08     0.66     0.98 1.00     2358     2494

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     5.05      1.34     2.80     8.02 1.00     2364     2146

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

• in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 19 raw observations.
• the estimated intercept (expected number of individuals for a mussel clump of log-transformed area = 0) is -1.1. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.33. When the mussel clump is 0, we expect to find 0.33 individuals in the mussel clump (which obviously does not exist if it has a log-transformed area of 0!).
  
• the estimated slope (rate at which the number of individuals changes per 1 unit change in log-transformed area), is 0.82 (mean) or 0.98 (median) with a standard deviation of 0.08. The 90% credibility intervals indicate that we are 90% confident that the slope is between 0.82 and 1 - e.g. there is a significant positive trend. When back-transformed onto the reponse scale, we see that for a one unit increase in log-transformed mussel clump area, the number of individuals increases by a factor of 2.27. This represents a 127% increase per 1 unit change in log-transformed mussel clump area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

summarise_draws (posterior)

peake.brm4 |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 9 × 8
  variable           median     lower    upper  rhat length ess_bulk ess_tail
  <chr>               <dbl>     <dbl>    <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept       -1.11     -2.44      0.278  1.00   2400    2401.    2461.
2 b_logAREA          0.819     0.662     0.983  1.00   2400    2358.    2494.
3 shape              4.92      2.64      7.74   1.00   2400    2364.    2146.
4 Intercept          5.73      5.56      5.92   1.00   2400    2430.    2344.
5 prior_Intercept    5.97      3.19      9.00   1.00   2400    2435.    2368.
6 prior_b           -0.0192   -1.87      2.01   1.00   2400    2467.    2277.
7 prior_shape        1.62      0.0490    4.71   1.00   2400    2445.    2258.
8 lprior            -5.93     -8.26     -4.22   1.00   2400    2371.    1895.
9 lp__            -159.     -162.     -158.     1.00   2400    2396.    2228.

We can also alter the CI level.

peake.brm4 |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x, credMass = 0.9),
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 9 × 8
  variable           median     lower     upper  rhat length ess_bulk ess_tail
  <chr>               <dbl>     <dbl>     <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept       -1.11     -2.18      0.0733  1.00   2400    2401.    2461.
2 b_logAREA          0.819     0.687     0.956   1.00   2400    2358.    2494.
3 shape              4.92      2.87      7.11    1.00   2400    2364.    2146.
4 Intercept          5.73      5.58      5.88    1.00   2400    2430.    2344.
5 prior_Intercept    5.97      3.41      8.25    1.00   2400    2435.    2368.
6 prior_b           -0.0192   -1.65      1.65    1.00   2400    2467.    2277.
7 prior_shape        1.62      0.0439    3.83    1.00   2400    2445.    2258.
8 lprior            -5.93     -7.64     -4.33    1.00   2400    2371.    1895.
9 lp__            -159.     -161.     -158.      1.00   2400    2396.    2228.

To narrow down to just the parameters of interest, see the code under the tidy_draws tab.

tidyMCMC

tidyMCMC(peake.brm4$fit, estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval", rhat = TRUE, ess = TRUE)

# A tibble: 8 × 7
  term            estimate std.error conf.low conf.high  rhat   ess
  <chr>              <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
1 b_Intercept      -1.10      0.695   -2.44       0.278 1.00   2388
2 b_logAREA         0.818     0.0821   0.662      0.983 1.00   2349
3 shape             5.05      1.34     2.64       7.74  0.999  2349
4 Intercept         5.73      0.0926   5.56       5.92  1.00   2424
5 prior_Intercept   5.96      1.49     3.19       9.00  1.000  2426
6 prior_b          -0.0141    1.01    -1.87       2.01  1.000  2472
7 prior_shape       1.96      1.41     0.0490     4.71  1.00   2263
8 lprior           -6.06      1.08    -8.26      -4.22  0.999  2349

Conclusions:

• the estimated intercept (expected number of individuals for a mussel clump of log-transformed area = 0) is -1.1. This is the median of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.33. When the mussel clump is 0, we expect to find 0.33 individuals in the mussel clump (which obviously does not exist if it has a log-transformed area of 0!).
  
• the estimated slope (rate at which the number of individuals changes per 1 unit change in log-transformed area), is 0.82 (mean) with a standard error of 0.08. The 95% credibility intervals indicate that we are 90% confident that the slope is between 0.66 and 0.98 - e.g. there is a significant positive trend. When back-transformed onto the reponse scale, we see that for a one unit increase in log-transformed mussel clump area, the number of individuals increases by a factor of 2.27. This represents a 127% increase per 1 unit change in log-transformed mussel clump area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

as_draws_df (posterior)

peake.brm4 |> as_draws_df()

# A draws_df: 800 iterations, 3 chains, and 9 variables
   b_Intercept b_logAREA shape Intercept prior_Intercept prior_b prior_shape
1        -0.35      0.74   4.4       5.9             9.8   -0.01        1.10
2        -0.55      0.76   5.7       5.8             6.3    0.16        2.24
3        -0.84      0.79   8.0       5.7             5.1    1.04        1.41
4        -1.11      0.81   3.9       5.7             6.4    2.61        5.05
5        -0.48      0.74   4.0       5.7             7.2    1.12        1.43
6         0.44      0.64   4.8       5.8             5.2    1.37        1.63
7        -0.40      0.74   2.0       5.8             6.7   -1.74        0.80
8        -1.92      0.92   4.6       5.7             5.2    0.29        2.98
9        -0.78      0.79   4.5       5.8             5.7   -0.90        0.72
10       -0.57      0.74   7.3       5.6             6.5    0.64        3.22
   lprior
1    -5.4
2    -6.5
3    -8.5
4    -5.1
5    -5.1
6    -5.7
7    -3.8
8    -5.7
9    -5.5
10   -7.9
# ... with 2390 more draws, and 1 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

## summarised
peake.brm4 |>
  as_draws_df() |>
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    length,
    rhat,
    ess_bulk, ess_tail
  )

# A tibble: 9 × 8
  variable           median     lower    upper length  rhat ess_bulk ess_tail
  <chr>               <dbl>     <dbl>    <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept       -1.11     -2.44      0.278   2400  1.00    2401.    2461.
2 b_logAREA          0.819     0.662     0.983   2400  1.00    2358.    2494.
3 shape              4.92      2.64      7.74    2400  1.00    2364.    2146.
4 Intercept          5.73      5.56      5.92    2400  1.00    2430.    2344.
5 prior_Intercept    5.97      3.19      9.00    2400  1.00    2435.    2368.
6 prior_b           -0.0192   -1.87      2.01    2400  1.00    2467.    2277.
7 prior_shape        1.62      0.0490    4.71    2400  1.00    2445.    2258.
8 lprior            -5.93     -8.26     -4.22    2400  1.00    2371.    1895.
9 lp__            -159.     -162.     -158.      2400  1.00    2396.    2228.

## summarised on fractional scale
peake.brm4 |>
  as_draws_df() |>
  dplyr::select(starts_with("b_")) |>
  ## mutate(across(everything(), exp)) |>
  exp() |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    length,
    ess_bulk, ess_tail
  )

Warning: Dropping 'draws_df' class as required metadata was removed.

# A tibble: 2 × 8
  variable    median  lower upper  rhat length ess_bulk ess_tail
  <chr>        <dbl>  <dbl> <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept  0.330 0.0322  1.05 1.000   2400    2384.    2448.
2 b_logAREA    2.27  1.90    2.63 1.00    2400    2338.    2431.

tidy_draws

Return the draws (samples) for each parameter in wide format

peake.brm4 |> tidy_draws()

# A tibble: 2,400 × 18
   .chain .iteration .draw b_Intercept b_logAREA shape Intercept prior_Intercept
    <int>      <int> <int>       <dbl>     <dbl> <dbl>     <dbl>           <dbl>
 1      1          1     1      -0.346     0.743  4.38      5.86            9.80
 2      1          2     2      -0.551     0.762  5.67      5.81            6.26
 3      1          3     3      -0.843     0.787  7.98      5.73            5.11
 4      1          4     4      -1.11      0.814  3.91      5.69            6.38
 5      1          5     5      -0.480     0.743  3.97      5.72            7.23
 6      1          6     6       0.439     0.643  4.75      5.81            5.21
 7      1          7     7      -0.403     0.740  1.99      5.78            6.68
 8      1          8     8      -1.92      0.917  4.56      5.73            5.17
 9      1          9     9      -0.776     0.792  4.48      5.84            5.73
10      1         10    10      -0.568     0.739  7.33      5.60            6.50
# ℹ 2,390 more rows
# ℹ 10 more variables: prior_b <dbl>, prior_shape <dbl>, lprior <dbl>,
#   lp__ <dbl>, accept_stat__ <dbl>, stepsize__ <dbl>, treedepth__ <dbl>,
#   n_leapfrog__ <dbl>, divergent__ <dbl>, energy__ <dbl>

peake.brm4 |> get_variables()

 [1] "b_Intercept"     "b_logAREA"       "shape"           "Intercept"      
 [5] "prior_Intercept" "prior_b"         "prior_shape"     "lprior"         
 [9] "lp__"            "accept_stat__"   "stepsize__"      "treedepth__"    
[13] "n_leapfrog__"    "divergent__"     "energy__"       

peake.brm4$fit |>
  tidy_draws() |>
  dplyr::select(matches("^b_.*")) |>
  exp() |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 2 × 8
  variable    median  lower upper  rhat length ess_bulk ess_tail
  <chr>        <dbl>  <dbl> <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept  0.330 0.0322  1.05 1.000   2400    2384.    2448.
2 b_logAREA    2.27  1.90    2.63 1.00    2400    2338.    2431.

The above can be useful for exploring the full posteriors of all the parameters of performing specific calculations using the posteriors, summarising the parameters takes a few more steps.

peake.brm4 |>
  tidy_draws() |>
  exp() |>
  gather_variables() |>
  median_hdci()

# A tibble: 15 × 7
   .variable         .value   .lower   .upper .width .point .interval
   <chr>              <dbl>    <dbl>    <dbl>  <dbl> <chr>  <chr>    
 1 accept_stat__   2.59e+ 0 1.97e+ 0 2.72e+ 0   0.95 median hdci     
 2 b_Intercept     3.30e- 1 5.40e- 2 1.07e+ 0   0.95 median hdci     
 3 b_logAREA       2.27e+ 0 1.91e+ 0 2.64e+ 0   0.95 median hdci     
 4 divergent__     1   e+ 0 1   e+ 0 1   e+ 0   0.95 median hdci     
 5 energy__        8.16e+69 6.42e+68 3.09e+71   0.95 median hdci     
 6 Intercept       3.09e+ 2 2.59e+ 2 3.71e+ 2   0.95 median hdci     
 7 lp__            5.72e-70 3.48e-73 1.57e-69   0.95 median hdci     
 8 lprior          2.66e- 3 2.15e- 5 1.08e- 2   0.95 median hdci     
 9 n_leapfrog__    2.01e+ 1 2.01e+ 1 1.10e+ 3   0.95 median hdci     
10 prior_b         9.81e- 1 2.31e- 2 5.24e+ 0   0.95 median hdci     
11 prior_Intercept 3.92e+ 2 2.75e+ 0 4.50e+ 3   0.95 median hdci     
12 prior_shape     5.08e+ 0 1.04e+ 0 1.10e+ 2   0.95 median hdci     
13 shape           1.37e+ 2 5.11e+ 0 1.71e+ 3   0.95 median hdci     
14 stepsize__      2.11e+ 0 1.88e+ 0 2.26e+ 0   0.95 median hdci     
15 treedepth__     7.39e+ 0 2.72e+ 0 2.01e+ 1   0.95 median hdci     

The gather_draws on the other hand, conveniently combines tidy_draws and gather_variables together in a single command. Furthermore, it returns all of the variables. The spread_draws function allows users to target specific variables (either by naming them in full or via regexp).

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

peake.brm4 |> get_variables()

 [1] "b_Intercept"     "b_logAREA"       "shape"           "Intercept"      
 [5] "prior_Intercept" "prior_b"         "prior_shape"     "lprior"         
 [9] "lp__"            "accept_stat__"   "stepsize__"      "treedepth__"    
[13] "n_leapfrog__"    "divergent__"     "energy__"       

peake.draw <- peake.brm4 |> gather_draws(`b_.*`, regex = TRUE)
peake.draw

# A tibble: 4,800 × 5
# Groups:   .variable [2]
   .chain .iteration .draw .variable   .value
    <int>      <int> <int> <chr>        <dbl>
 1      1          1     1 b_Intercept -0.346
 2      1          2     2 b_Intercept -0.551
 3      1          3     3 b_Intercept -0.843
 4      1          4     4 b_Intercept -1.11 
 5      1          5     5 b_Intercept -0.480
 6      1          6     6 b_Intercept  0.439
 7      1          7     7 b_Intercept -0.403
 8      1          8     8 b_Intercept -1.92 
 9      1          9     9 b_Intercept -0.776
10      1         10    10 b_Intercept -0.568
# ℹ 4,790 more rows

We can then summarise this

peake.draw |> median_hdci()

# A tibble: 2 × 7
  .variable   .value .lower .upper .width .point .interval
  <chr>        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_Intercept -1.11  -2.44   0.287   0.95 median hdci     
2 b_logAREA    0.819  0.662  0.983   0.95 median hdci     

peake.brm4 |>
  gather_draws(`b_.*`, regex = TRUE) |>
  mutate(.value = exp(.value)) |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat, length, ess_bulk, ess_tail
  )

# A tibble: 2 × 9
# Groups:   .variable [2]
  .variable   variable median  lower upper  rhat length ess_bulk ess_tail
  <chr>       <chr>     <dbl>  <dbl> <dbl> <dbl>  <dbl>    <dbl>    <dbl>
1 b_Intercept .value    0.330 0.0322  1.05  1.00   2400    2401.    2461.
2 b_logAREA   .value    2.27  1.90    2.63  1.00   2400    2358.    2494.

peake.brm4 |>
  gather_draws(`b_.*`, regex = TRUE) |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

We could alternatively express the parameters on the response scale.

peake.brm4 |>
  gather_draws(`.Intercept.*|.*AREA.*`, regex = TRUE) |>
  group_by(.variable) |>
  mutate(.value = exp(.value)) |>
  median_hdci()

# A tibble: 1 × 7
  .variable .value .lower .upper .width .point .interval
  <chr>      <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_logAREA   2.27   1.91   2.64   0.95 median hdci     

Conclusions:

• the estimated intercept (expected number of individuals for a mussel clump of log-transformed area = 0) is -1.11. This is the median of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.33. When the mussel clump is 0, we expect to find 0.33 individuals in the mussel clump (which obviously does not exist if it has a log-transformed area of 0!).
  
• the estimated slope (rate at which the number of individuals changes per 1 unit change in log-transformed area), is 0.82 (mean) with a standard error of 0.66. The 95% credibility intervals indicate that we are 90% confident that the slope is between 0.98 and 0.95 - e.g. there is a significant positive trend. When back-transformed onto the reponse scale, we see that for a one unit increase in log-transformed mussel clump area, the number of individuals increases by a factor of 2.27. This represents a 127% increase per 1 unit change in log-transformed mussel clump area.
• Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

bayesplot

peake.brm4 |> plot(plotfun = "mcmc_intervals")

spread_draws

peake.brm4 |> spread_draws(`b_.*`, regex = TRUE)

# A tibble: 2,400 × 5
   .chain .iteration .draw b_Intercept b_logAREA
    <int>      <int> <int>       <dbl>     <dbl>
 1      1          1     1      -0.346     0.743
 2      1          2     2      -0.551     0.762
 3      1          3     3      -0.843     0.787
 4      1          4     4      -1.11      0.814
 5      1          5     5      -0.480     0.743
 6      1          6     6       0.439     0.643
 7      1          7     7      -0.403     0.740
 8      1          8     8      -1.92      0.917
 9      1          9     9      -0.776     0.792
10      1         10    10      -0.568     0.739
# ℹ 2,390 more rows

## summarised
peake.brm4 |>
  as_draws_df() |>
  dplyr::select(starts_with("b_")) |>
  summarise_draws(
    "median",
    ~ HDInterval::hdi(.x),
    "rhat",
    "ess_bulk"
  )

Warning: Dropping 'draws_df' class as required metadata was removed.

# A tibble: 2 × 6
  variable    median  lower upper  rhat ess_bulk
  <chr>        <dbl>  <dbl> <dbl> <dbl>    <dbl>
1 b_Intercept -1.11  -2.44  0.278 1.000    2384.
2 b_logAREA    0.819  0.662 0.983 1.00     2338.

## summarised on fractional scale
peake.brm4 |>
  as_draws_df() |>
  dplyr::select(starts_with("b_")) |>
  mutate(across(everything(), exp)) |>
  summarise_draws(
    "median",
    ~ HDInterval::hdi(.x),
    "rhat",
    "ess_bulk"
  )

Warning: Dropping 'draws_df' class as required metadata was removed.

# A tibble: 2 × 6
  variable    median  lower upper  rhat ess_bulk
  <chr>        <dbl>  <dbl> <dbl> <dbl>    <dbl>
1 b_Intercept  0.330 0.0322  1.05 1.000    2384.
2 b_logAREA    2.27  1.90    2.63 1.00     2338.

posterior_samples

peake.brm4 |>
  posterior_samples() |>
  as_tibble()

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

# A tibble: 2,400 × 9
   b_Intercept b_logAREA shape Intercept prior_Intercept  prior_b prior_shape
         <dbl>     <dbl> <dbl>     <dbl>           <dbl>    <dbl>       <dbl>
 1      -0.346     0.743  4.38      5.86            9.80 -0.01000       1.10 
 2      -0.551     0.762  5.67      5.81            6.26  0.163         2.24 
 3      -0.843     0.787  7.98      5.73            5.11  1.04          1.41 
 4      -1.11      0.814  3.91      5.69            6.38  2.61          5.05 
 5      -0.480     0.743  3.97      5.72            7.23  1.12          1.43 
 6       0.439     0.643  4.75      5.81            5.21  1.37          1.63 
 7      -0.403     0.740  1.99      5.78            6.68 -1.74          0.796
 8      -1.92      0.917  4.56      5.73            5.17  0.294         2.98 
 9      -0.776     0.792  4.48      5.84            5.73 -0.902         0.716
10      -0.568     0.739  7.33      5.60            6.50  0.644         3.22 
# ℹ 2,390 more rows
# ℹ 2 more variables: lprior <dbl>, lp__ <dbl>

\(R^2\)

peake.brm4 |>
  bayes_R2(summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.7334567 0.6723578 0.7518716   0.95 median      hdci

Modelsummary

peake.brm4 |> modelsummary(
  statistic = c("conf.low", "conf.high"),
  shape = term ~ statistic
)

Warning: 
`modelsummary` uses the `performance` package to extract goodness-of-fit
statistics from models of this class. You can specify the statistics you wish
to compute by supplying a `metrics` argument to `modelsummary`, which will then
push it forward to `performance`. Acceptable values are: "all", "common",
"none", or a character vector of metrics names. For example: `modelsummary(mod,
metrics = c("RMSE", "R2")` Note that some metrics are computationally
expensive. See `?performance::performance` for details.
 This warning appears once per session.

	(1)
	Est.	2.5 %	97.5 %
b_Intercept	-1.108	-2.470	0.276
b_logAREA	0.819	0.656	0.979
Num.Obs.	25
R2	0.733
ELPD	-156.4
ELPD s.e.	5.4
LOOIC	312.8
LOOIC s.e.	10.8
WAIC	312.7
RMSE	249.40

peake.brm4 |> modelsummary(
  statistic = c("conf.low", "conf.high"),
  shape = term ~ statistic,
  exponentiate = TRUE
)

	(1)
	Est.	2.5 %	97.5 %
b_Intercept	0.330	0.085	1.318
b_logAREA	2.269	1.927	2.662
Num.Obs.	25
R2	0.733
ELPD	-156.4
ELPD s.e.	5.4
LOOIC	312.8
LOOIC s.e.	10.8
WAIC	312.7
RMSE	246.41

modelsummary(list("Raw" = peake.brm4, "Exponentiated" = peake.brm4),
  statistic = c("conf.low", "conf.high"),
  shape = term ~ statistic,
  exponentiate = c(FALSE, TRUE)
)

	Raw			Exponentiated
	Est.	2.5 %	97.5 %	Est.	2.5 %	97.5 %
b_Intercept	-1.108	-2.470	0.276	0.330	0.085	1.318
b_logAREA	0.819	0.656	0.979	2.269	1.927	2.662
Num.Obs.	25			25
R2	0.733			0.733
ELPD	-156.4			-156.4
ELPD s.e.	5.4			5.4
LOOIC	312.8			312.8
LOOIC s.e.	10.8			10.8
WAIC	312.7			312.7
RMSE	250.52			243.99

peake.brm4 |> modelplot()

peake.brm4 |> modelplot(exponentiate = TRUE)

11 Hypothesis testing

rstanarm

Prob. of effect

Unfortunately, the inline log-transformation of the predictor interfers with hypothesis()’s ability to find the slope. We can get around this by renaming before calling hypothesis().

peake.rstanarm4 |>
  as.data.frame() |>
  rename(lAREA = `log(AREA)`) |>
  hypothesis("lAREA>0")

Hypothesis Tests for class :
   Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (lAREA) > 0     0.82      0.08     0.68     0.96        Inf         1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Conclusions:

• the probability of a positive effect of mussel clump area on number of individuals is 1. That is, there is, we are 100% confident that there is an effect.
• the evidence ratio is normally the ratio of the number of cases that satisty the hypothesis to the number of cases that do not. Since the number of cases that do not satisfy the hypothesis is 0, the evidence ratio is Inf (since division by 0)

Alternatively…

peake.rstanarm4 |>
  tidy_draws() |>
  summarise(P = mean(`log(AREA)` > 0))

# A tibble: 1 × 1
      P
  <dbl>
1     1

Conclusions:

• the probability of a positive effect of mussel clump area on number of individuals is 1. That is, there is, we are 100% confident that there is an effect.

Change 5000 to 10000 (Method 1)

newdata <- list(AREA = c(5000, 10000))
## fractional scale
peake.rstanarm4 |>
  emmeans(~AREA, at = newdata, type = "response") |>
  pairs(reverse = TRUE)

 contrast             ratio lower.HPD upper.HPD
 AREA10000 / AREA5000  1.77      1.56      1.96

Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

## absolute scale
peake.rstanarm4 |>
  emmeans(~AREA, at = newdata) |>
  regrid() |>
  pairs(reverse = TRUE)

 contrast             estimate lower.HPD upper.HPD
 AREA10000 - AREA5000      272       184       369

Point estimate displayed: median 
HPD interval probability: 0.95 

Conclusions:

• the change in expected number of individuals associated with an increase in mussel clump area from 5,000 to 10,000 is 0.57.

If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.

peake.mcmc <- peake.rstanarm4 |>
  emmeans(~AREA, at = newdata) |>
  regrid() |>
  tidy_draws() |>
  rename_with(~ str_replace(., "AREA ", "p")) |>
  mutate(
    Eff = p10000 - p5000,
    PEff = 100 * Eff / p5000
  )
peake.mcmc |> head()

# A tibble: 6 × 7
  .chain .iteration .draw p5000 p10000   Eff  PEff
   <int>      <int> <int> <dbl>  <dbl> <dbl> <dbl>
1      1          1     1  350.   569.  220.  62.8
2      1          2     2  363.   605.  242.  66.5
3      1          3     3  395.   695.  300.  76.0
4      1          4     4  353.   567.  214.  60.5
5      1          5     5  299.   607.  309. 103. 
6      1          6     6  418.   767.  349.  83.5

Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).

peake.mcmc |> tidyMCMC(
  estimate.method = "median",
  conf.int = TRUE, conf.method = "HPDinterval"
)

# A tibble: 7 × 5
  term       estimate std.error conf.low conf.high
  <chr>         <dbl>     <dbl>    <dbl>     <dbl>
1 .chain          2       0.817      1         3  
2 .iteration    400.    231.         1       761  
3 .draw        1200.    693.         1      2281  
4 p5000         356.     35.6      290.      427. 
5 p10000        632.     76.4      495.      782. 
6 Eff           276.     48.8      184.      369. 
7 PEff           77.3    10.5       55.7      96.3

Alternatively, we could use median_hdci

peake.mcmc |> median_hdci(PEff)

# A tibble: 1 × 6
   PEff .lower .upper .width .point .interval
  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1  77.1   56.4   97.3   0.95 median hdci     

Conclusions:

• the percentage change in expected number of individuals associated with an increase in mussel clump area from 5,000 to 10,000 is 77.07% with a 95% credibility interval of 56.43 97.26

To get the probability that the effect is greater than a 50% increase.

peake.mcmc |> summarise(P = mean(PEff > 50))

# A tibble: 1 × 1
      P
  <dbl>
1 0.998

Conclusions:

• the probability that the number of individuals will increase by more than 50% following an increase in mussel clump area from 5,000 to 10,000 is

Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).

peake.mcmc |> hypothesis("PEff>50")

Hypothesis Tests for class :
       Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (PEff)-(50) > 0    27.31     10.47    10.51     44.7        479         1
  Star
1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Change 5000 to 10000 (Method 2)

## Note, the P value is on absolute difference
newdata <- list(AREA = c(5000, 10000))
peake.rstanarm4 |>
  emmeans(~AREA, at = newdata) |>
  regrid() |>
  pairs(reverse = TRUE) |>
  tidy_draws() |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)

# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1     1  272.  184.  369.   0.95 median hdci     

We can also express this as a percentage change

newdata <- list(AREA = c(5000, 10000))
## Simple
peake.rstanarm4 |>
  emmeans(~AREA, at = newdata) |>
  pairs(reverse = TRUE) |>
  regrid()

 contrast           ratio lower.HPD upper.HPD
 AREA10000/AREA5000  1.77      1.56      1.96

Point estimate displayed: median 
HPD interval probability: 0.95 

## More advanced (both P and percent change)
peake.mcmc <- peake.rstanarm4 |>
  emmeans(~AREA, at = newdata) |>
  pairs(reverse = TRUE) |>
  regrid() |>
  tidy_draws() |>
  mutate(across(contains("contrast"), ~ 100 * (. - 1)))

peake.mcmc |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)

# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 0.998  77.1  56.4  97.3   0.95 median hdci     

Change 5000 to 10000 (Method 3)

peake.rstanarm4 |>
  linpred_draws(newdata = as.data.frame(newdata)) |>
  mutate(.linpred = exp(.linpred)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = diff(.linpred),
    PEff = 100 * Eff / .linpred[1]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()

# A tibble: 3 × 7
  name    value  .lower  .upper .width .point .interval
  <chr>   <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   272.    184.    369.      0.95 median hdci     
2 P       0.998   0.998   0.998   0.95 median hdci     
3 PEff   77.1    56.4    97.3     0.95 median hdci     

## OR
peake.rstanarm4 |>
  epred_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = diff(.epred),
    PEff = 100 * Eff / .epred[1]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()

# A tibble: 3 × 7
  name    value  .lower  .upper .width .point .interval
  <chr>   <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   272.    184.    369.      0.95 median hdci     
2 P       0.998   0.998   0.998   0.95 median hdci     
3 PEff   77.1    56.4    97.3     0.95 median hdci     

## OR for prediction of individual values
peake.rstanarm4 |>
  predicted_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = diff(.prediction),
    PEff = 100 * Eff / .prediction[1]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()

# A tibble: 3 × 7
  name    value   .lower   .upper .width .point .interval
  <chr>   <dbl>    <dbl>    <dbl>  <dbl> <chr>  <chr>    
1 Eff   234     -371.    1018.      0.95 median hdci     
2 P       0.590    0.590    0.590   0.95 median hdci     
3 PEff   75.2    -82.3    495.      0.95 median hdci     

brm

Prob. of effect

peake.brm4 |> hypothesis("logAREA>0")

Hypothesis Tests for class b:
     Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (logAREA) > 0     0.82      0.08     0.68     0.95        Inf         1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Conclusions:

• the probability of a positive effect of mussel clump area on number of individuals is 1. That is, there is, we are 100% confident that there is an effect.
• the evidence ratio is normally the ratio of the number of cases that satisty the hypothesis to the number of cases that do not. Since the number of cases that do not satisfy the hypothesis is 0, the evidence ratio is Inf (since division by 0)

Alternatively…

peake.brm4 |>
  tidy_draws() |>
  summarise(P = mean(b_logAREA > 0))

# A tibble: 1 × 1
      P
  <dbl>
1     1

Conclusions:

• the probability of a positive effect of mussel clump area on number of individuals is 1. That is, there is, we are 100% confident that there is an effect.

Change 5000 to 10000 (Method 1)

newdata <- list(AREA = c(5000, 10000))
## fractional scale
peake.brm4 |>
  emmeans(~AREA, at = newdata, type = "response") |>
  pairs(reverse = TRUE)

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

 contrast             ratio lower.HPD upper.HPD
 AREA10000 / AREA5000  1.76      1.56      1.96

Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

## absolute scale
peake.brm4 |>
  emmeans(~AREA, at = newdata) |>
  regrid() |>
  pairs(reverse = TRUE)

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

 contrast             estimate lower.HPD upper.HPD
 AREA10000 - AREA5000      268       189       367

Point estimate displayed: median 
HPD interval probability: 0.95 

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

Conclusions:

• the change in expected number of individuals associated with an increase in mussel clump area from 5,000 to 10,000 is 0.57.

If we want to derive other properties, such as the percentage change, then we use tidy_draws() and then simple tidyverse spreadsheet operation.

peake.mcmc <- peake.brm4 |>
  emmeans(~AREA, at = newdata) |>
  regrid() |>
  tidy_draws() |>
  rename_with(~ str_replace(., "AREA ", "p")) |>
  mutate(
    Eff = p10000 - p5000,
    PEff = 100 * Eff / p5000
  )

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

peake.mcmc |> head()

# A tibble: 6 × 7
  .chain .iteration .draw p5000 p10000   Eff  PEff
   <int>      <int> <int> <dbl>  <dbl> <dbl> <dbl>
1      1          1     1  395.   661.  266.  67.3
2      1          2     2  379.   643.  264.  69.6
3      1          3     3  350.   604.  254.  72.5
4      1          4     4  337.   592.  255.  75.8
5      1          5     5  345.   578.  232.  67.3
6      1          6     6  372.   581.  209.  56.2

Now we can calculate the medians and HPD intervals of each column (and ignore the .chain, .iteration and .draw).

peake.mcmc |> tidyMCMC(
  estimate.method = "median",
  conf.int = TRUE, conf.method = "HPDinterval"
)

# A tibble: 7 × 5
  term       estimate std.error conf.low conf.high
  <chr>         <dbl>     <dbl>    <dbl>     <dbl>
1 .chain          2       0.817      1         3  
2 .iteration    400.    231.         1       761  
3 .draw        1200.    693.         1      2281  
4 p5000         355.     33.1      290.      418. 
5 p10000        627.     72.4      491.      768. 
6 Eff           272.     46.6      189.      367. 
7 PEff           76.6    10.1       56.2      95.7

Alternatively, we could use median_hdci

peake.mcmc |> median_hdci(PEff)

# A tibble: 1 × 6
   PEff .lower .upper .width .point .interval
  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1  76.4   56.3   95.8   0.95 median hdci     

Conclusions:

• the percentage change in expected number of individuals associated with an increase in mussel clump area from 5,000 to 10,000 is 76.44% with a 95% credibility interval of 56.29 95.81

To get the probability that the effect is greater than a 50% increase.

peake.mcmc |> summarise(P = mean(PEff > 50))

# A tibble: 1 × 1
      P
  <dbl>
1 0.997

Conclusions:

• the probability that the number of individuals will increase by more than 50% following an increase in mussel clump area from 5,000 to 10,000 is

Finally, we could alternatively use hypothesis(). Note that when we do so, the estimate is the difference between the effect and the hypothesised effect (50%).

peake.mcmc |> hypothesis("PEff>50")

Hypothesis Tests for class :
       Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (PEff)-(50) > 0    26.61     10.07     10.6    43.64     341.86         1
  Star
1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Change 5000 to 10000 (Method 2)

## Note, the P value is on absolute difference
newdata <- list(AREA = c(5000, 10000))
peake.brm4 |>
  emmeans(~AREA, at = newdata) |>
  regrid() |>
  pairs(reverse = TRUE) |>
  tidy_draws() |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1     1  268.  189.  367.   0.95 median hdci     

We can also express this as a percentage change

newdata <- list(AREA = c(5000, 10000))
## Simple
peake.brm4 |>
  emmeans(~AREA, at = newdata) |>
  pairs(reverse = TRUE) |>
  regrid()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

 contrast           ratio lower.HPD upper.HPD
 AREA10000/AREA5000  1.76      1.56      1.96

Point estimate displayed: median 
HPD interval probability: 0.95 

## More advanced (both P and percent change)
peake.mcmc <- peake.brm4 |>
  emmeans(~AREA, at = newdata) |>
  pairs(reverse = TRUE) |>
  regrid() |>
  tidy_draws() |>
  mutate(across(contains("contrast"), ~ 100 * (. - 1)))

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

peake.mcmc |>
  summarise(across(contains("contrast"),
    list(
      P = ~ mean(. > 50),
      HDCI = ~ median_hdci(.)
    ),
    .names = c("{.fn}")
  )) |>
  tidyr::unpack(HDCI)

# A tibble: 1 × 7
      P     y  ymin  ymax .width .point .interval
  <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 0.997  76.4  56.3  95.8   0.95 median hdci     

Change 5000 to 10000 (Method 3)

peake.brm4 |>
  linpred_draws(newdata = as.data.frame(newdata)) |>
  mutate(.linpred = exp(.linpred)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = diff(.linpred),
    PEff = 100 * Eff / .linpred[1]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()

# A tibble: 3 × 7
  name    value  .lower  .upper .width .point .interval
  <chr>   <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   268.    189.    367.      0.95 median hdci     
2 P       0.997   0.997   0.997   0.95 median hdci     
3 PEff   76.4    56.3    95.8     0.95 median hdci     

## OR
peake.brm4 |>
  epred_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = diff(.epred),
    PEff = 100 * Eff / .epred[1]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()

# A tibble: 3 × 7
  name    value  .lower  .upper .width .point .interval
  <chr>   <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 Eff   268.    189.    367.      0.95 median hdci     
2 P       0.997   0.997   0.997   0.95 median hdci     
3 PEff   76.4    56.3    95.8     0.95 median hdci     

## OR for prediction of individual values
peake.brm4 |>
  predicted_draws(newdata = as.data.frame(newdata)) |>
  ungroup() |>
  group_by(.draw) |>
  summarise(
    Eff = diff(.prediction),
    PEff = 100 * Eff / .prediction[1]
  ) |>
  ungroup() |>
  mutate(P = mean(PEff > 50)) |>
  pivot_longer(cols = -.draw) |>
  group_by(name) |>
  median_hdci()

# A tibble: 3 × 7
  name    value   .lower   .upper .width .point .interval
  <chr>   <dbl>    <dbl>    <dbl>  <dbl> <chr>  <chr>    
1 Eff   236.    -371.    1007.      0.95 median hdci     
2 P       0.598    0.598    0.598   0.95 median hdci     
3 PEff   74.7    -81.2    452.      0.95 median hdci     

12 Summary figure

rstanarm

## Using emmeans
peake.grid <- with(peake, list(AREA = seq(min(AREA), max(AREA), len = 100)))

newdata <- emmeans(peake.rstanarm4, ~AREA, at = peake.grid, type = "response") |> as.data.frame()
head(newdata)

      AREA      prob lower.HPD upper.HPD
  462.2500  49.39540  31.75184  74.76083
  731.7629  72.31525  49.50751 102.29683
 1001.2759  93.80586  65.63684 125.19252
 1270.7888 114.38392  83.50520 148.26200
 1540.3017 133.95323  99.45183 170.43119
 1809.8146 152.95173 117.95764 193.16183

Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

ggplot(newdata, aes(y = prob, x = AREA)) +
  geom_point(data = peake, aes(y = INDIV)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
  scale_y_continuous("Individuals") +
  scale_x_continuous("Mussel clump area") +
  theme_classic()

## spaghetti plot
newdata <- emmeans(peake.rstanarm4, ~AREA, at = peake.grid) |>
  regrid() |>
  gather_emmeans_draws()
newdata |> head()

# A tibble: 6 × 5
# Groups:   AREA [1]
   AREA .chain .iteration .draw .value
  <dbl>  <int>      <int> <int>  <dbl>
1  462.     NA         NA     1   65.6
2  462.     NA         NA     2   63.1
3  462.     NA         NA     3   56.6
4  462.     NA         NA     4   69.6
5  462.     NA         NA     5   26.0
6  462.     NA         NA     6   51.9

ggplot(newdata, aes(y = .value, x = AREA)) +
  geom_line(aes(group = .draw), colour = "orange", alpha = 0.01) +
  geom_point(data = peake, aes(y = INDIV)) +
  scale_y_continuous("Number of Individuals") +
  scale_x_continuous(expression(Mussel ~ clump ~ area ~ (per ~ mm^2))) +
  theme_classic()

brms

## Using emmeans
peake.grid <- with(peake, list(AREA = seq(min(AREA), max(AREA), len = 100)))

newdata <- emmeans(peake.brm4, ~AREA, at = peake.grid, type = "response") |> as.data.frame()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

head(newdata)

      AREA      prob lower.HPD upper.HPD
  462.2500  50.21763  31.10871  71.78707
  731.7629  73.02912  50.88985 100.65473
 1001.2759  94.37914  68.19245 124.33004
 1270.7888 114.99564  86.42094 147.87021
 1540.3017 134.72080 104.19582 169.97697
 1809.8146 153.53609 120.81942 190.86984

Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

ggplot(newdata, aes(y = prob, x = AREA)) +
  geom_point(data = peake, aes(y = INDIV)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD), fill = "blue", alpha = 0.3) +
  scale_y_continuous("Individuals") +
  scale_x_continuous("Mussel clump area") +
  theme_classic()

## spaghetti plot
newdata <- emmeans(peake.brm4, ~AREA, at = peake.grid) |>
  regrid() |>
  gather_emmeans_draws()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

newdata |> head()

# A tibble: 6 × 5
# Groups:   AREA [1]
   AREA .chain .iteration .draw .value
  <dbl>  <int>      <int> <int>  <dbl>
1  462.     NA         NA     1   67.4
2  462.     NA         NA     2   61.8
3  462.     NA         NA     3   53.8
4  462.     NA         NA     4   48.5
5  462.     NA         NA     5   58.9
6  462.     NA         NA     6   80.4

ggplot(newdata, aes(y = .value, x = AREA)) +
  geom_line(aes(group = .draw), colour = "orange", alpha = 0.01) +
  geom_point(data = peake, aes(y = INDIV)) +
  scale_y_continuous("Number of Individuals") +
  scale_x_continuous(expression(Mussel ~ clump ~ area ~ (per ~ mm^2))) +
  theme_classic()

13 References

Peake, A. J., and G. P. Quinn. 1993. “Temporal Variation in Species-Area Curves for Invertebrates in Clumps of an Intertidal Mussel.” Ecography 16: 269–77.